// Copyright 2022 Jeroen van Nugteren

// Permission is hereby granted, free of charge, to any person obtaining a copy of this software
// and associated documentation files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included in all copies or
// substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
// OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

// include header
#include "dfmultipole.hh"

// general headers
#include <functional>

// common headers
#include "rat/common/extra.hh"
#include "rat/common/marchingcubes.hh"

// code specific to Rat
namespace rat{namespace dm{


	// constructor
	DFHarmonic::DFHarmonic(){
		set_name("Harmonic");
	}

	DFHarmonic::DFHarmonic(
		const fltp reference_radius,
		const arma::uword num_poles, 
		const fltp flux_density, 
		const bool is_skew) : DFHarmonic(){
		set_reference_radius(reference_radius);
		set_num_poles(num_poles);
		set_flux_density(flux_density);
		set_is_skew(is_skew);
	}

	// factory
	ShDFHarmonicPr DFHarmonic::create(){
		return std::make_shared<DFHarmonic>();
	}

	ShDFHarmonicPr DFHarmonic::create(
		const fltp reference_radius,
		const arma::uword num_poles, 
		const fltp flux_density, 
		const bool is_skew){
		return std::make_shared<DFHarmonic>(reference_radius,num_poles,flux_density,is_skew);
	}

	// setters
	void DFHarmonic::set_num_poles(const arma::uword num_poles){
		num_poles_ = num_poles;
	}

	void DFHarmonic::set_flux_density(const fltp flux_density){
		flux_density_ = flux_density;
	}

	void DFHarmonic::set_is_skew(const bool is_skew){
		is_skew_ = is_skew;
	}

	void DFHarmonic::set_reference_radius(const fltp reference_radius){
		reference_radius_ = reference_radius;
	}



	// getters
	arma::uword DFHarmonic::get_num_poles()const{
		return num_poles_;
	}

	fltp DFHarmonic::get_flux_density()const{
		return flux_density_;
	}

	bool DFHarmonic::get_is_skew()const{
		return is_skew_;
	}

	fltp DFHarmonic::get_reference_radius()const{
		return reference_radius_;
	}

	// calculate coefficient
	std::complex<fltp> DFHarmonic::calc_coefficient()const{
		// create coefficient
		const fltp fref = (num_poles_+1)*std::pow(reference_radius_,num_poles_);
		return is_skew_ ? -(std::complex<fltp>(0.0,1.0)*flux_density_)/fref : -flux_density_/fref;
	}

	// validity check
	bool DFHarmonic::is_valid(const bool enable_throws)const{
		return true;
	}

	// type string for serialization
	std::string DFHarmonic::get_type(){
		return "rat::dm::dfharmonic";
	}

	// method for serialization into json
	void DFHarmonic::serialize(Json::Value &js, cmn::SList &list) const{
		// parent
		Node::serialize(js,list);

		// store type
		js["type"] = get_type();

	}

	// method for deserialisation from json
	void DFHarmonic::deserialize(
		const Json::Value &js, cmn::DSList &list, 
		const cmn::NodeFactoryMap &factory_list, 
		const boost::filesystem::path &pth){
		
		// parent
		Node::deserialize(js,list,factory_list,pth);
	}



	// constructcor
	DFMultipole::DFMultipole(){

	}

	DFMultipole::DFMultipole(
		const std::list<ShDFHarmonicPr>& harmonics){
		for(const auto& harmonic : harmonics)
			add_harmonic(harmonic);
	}

	// factory
	ShDFMultipolePr DFMultipole::create(){
		return std::make_shared<DFMultipole>();
	}

	ShDFMultipolePr DFMultipole::create(
		const std::list<ShDFHarmonicPr>& harmonics){
		return std::make_shared<DFMultipole>(harmonics);
	}


	// add a harmonic
	void DFMultipole::add_harmonic(const ShDFHarmonicPr& harmonic){
		harmonics_[harmonics_.size()+1] = harmonic;
	}

	// check if we have a pure dipole
	bool DFMultipole::is_pure_dipole()const{
		for(auto it=harmonics_.begin();it!=harmonics_.end();it++)
			if(it->second->get_num_poles()!=1)return false;
		return true;
	}

	// calculate field
	arma::Mat<std::complex<fltp> > DFMultipole::calc_F(
		const arma::Mat<std::complex<fltp> > &z)const{

		// calculate values
		arma::Mat<std::complex<fltp> > F(z.n_rows,z.n_cols,arma::fill::zeros);

		// walk over harmonics and add contributions
		for(auto it=harmonics_.begin();it!=harmonics_.end();it++){
			// get harmonic
			const ShDFHarmonicPr& harmonic = it->second;

			// calculate F
			const arma::uword np = harmonic->get_num_poles();
			F += harmonic->calc_coefficient()*arma::pow(z,np);
		}

		// return values
		return F; 
	}

	// calculate field
	std::complex<fltp> DFMultipole::calc_F(
		const std::complex<fltp> &z)const{

		// calculate values
		std::complex<fltp> F(RAT_CONST(0.0),RAT_CONST(0.0));

		// walk over harmonics and add contributions
		for(auto it=harmonics_.begin();it!=harmonics_.end();it++){
			// get harmonic
			const ShDFHarmonicPr& harmonic = it->second;

			// calculate F
			const arma::uword np = harmonic->get_num_poles();
			F += harmonic->calc_coefficient()*std::pow(z,np);
		}

		// return values
		return F; 
	}

	// calculate first derivative towards z
	arma::Mat<std::complex<fltp> > DFMultipole::calc_dFdz(
		const arma::Mat<std::complex<fltp> > &z)const{
		
		// calculate values
		arma::Mat<std::complex<fltp> > dF(z.n_rows,z.n_cols,arma::fill::zeros);

		// walk over harmonics and add contributions
		for(auto it=harmonics_.begin();it!=harmonics_.end();it++){
			// get harmonic
			const ShDFHarmonicPr& harmonic = it->second;
			const arma::uword np = harmonic->get_num_poles();
			dF += fltp(np+1)*harmonic->calc_coefficient()*arma::pow(z,arma::sword(np)-1);
		}

		// return values
		return dF; 
	}

	// calculate second derivative towards z
	arma::Mat<std::complex<fltp> > DFMultipole::calc_d2Fdz2(
		const arma::Mat<std::complex<fltp> > &z)const{

		// calculate values
		arma::Mat<std::complex<fltp> > ddF(z.n_rows,z.n_cols,arma::fill::zeros);

		// walk over harmonics and add contributions
		for(auto it=harmonics_.begin();it!=harmonics_.end();it++){
			// get harmonic
			const ShDFHarmonicPr& harmonic = it->second;
			const arma::uword np = harmonic->get_num_poles();
			ddF += fltp(np)*fltp(np+1)*harmonic->calc_coefficient()*arma::pow(z,arma::sword(np)-2);
		}

		// return values
		return ddF; 
	}

	// derivative of A to z
	arma::Mat<std::complex<fltp> > DFMultipole::calc_dAdz(
		const arma::Mat<std::complex<fltp> >&z)const{
		return arma::conj(calc_dFdz(z));
	}

	// derivative of V to z
	arma::Mat<std::complex<fltp> > DFMultipole::calc_dVdz(
		const arma::Mat<std::complex<fltp> >&z)const{
		return std::complex<fltp>(0.0,1.0)*calc_dAdz(z);
	}

	// flux density vector
	arma::Mat<std::complex<fltp> > DFMultipole::calc_B(
		const arma::Mat<std::complex<fltp> >&z)const{
		return -calc_dVdz(z);
	}

	// calculate flux density magnitude
	arma::Mat<fltp> DFMultipole::calc_Bmag(
		const arma::Mat<std::complex<fltp> >&z)const{
		return arma::abs(calc_B(z));
	}


	// create aperture in complex plane
	arma::Row<std::complex<fltp> > DFMultipole::create_aperture(
		const arma::Row<fltp>&t)const{
		// create theta
		const arma::Row<fltp> theta = 2*arma::Datum<fltp>::pi*t;

		// calculate coordinates
		const arma::Mat<std::complex<fltp> > Ra(
			aperture_a_*arma::pow(arma::abs(arma::cos(theta)),2.0/aperture_n_)%arma::sign(arma::cos(theta)),
			aperture_b_*arma::pow(arma::abs(arma::sin(theta)),2.0/aperture_n_)%arma::sign(arma::sin(theta)));

		// return aperture coordinates
		return Ra;
	}

	// create aperture in complex plane
	std::complex<fltp> DFMultipole::create_aperture(const fltp t)const{
		// create theta
		const fltp theta = 2*arma::Datum<fltp>::pi*t;

		// calculate coordinates
		const std::complex<fltp> Ra(
			aperture_a_*std::pow(std::abs(std::cos(theta)),2.0/aperture_n_)*cmn::Extra::sign(std::cos(theta)),
			aperture_b_*std::pow(std::abs(std::sin(theta)),2.0/aperture_n_)*cmn::Extra::sign(std::sin(theta)));

		// return aperture coordinates
		return Ra;
	}

	// create poles
	ShDFPolygonPr DFMultipole::create_poles()const{
		// find where the poles are tangent around the aperture
		const arma::Row<fltp> t = arma::linspace<arma::Row<fltp> >(RAT_CONST(0.0),RAT_CONST(1.0),1001llu);
		const arma::Row<std::complex<fltp> > z_aperture = create_aperture(t);


		std::cout<<calc_dAdz(z_aperture)<<std::endl;



		// // Ampere-turns per pole
		// const arma::Row<fltp> NI = RAT_CONST(1.0)/arma::Datum<fltp>::mu_0*arma::imag(calc_F(z_aperture)); // [A]

		// // find where the aperture is parallel to F
		// std::list<std::pair<fltp, fltp> > min_max;

		// // walk over time steps and search for local minima and maxima
		// for(arma::uword i=1;i<t.n_elem;i++){
		// 	if((NI(i)-NI(i-1))*(NI(i)-NI(i+1))>0)
		// 		min_max.push_back({t(i),NI(i)});
		// 	else if(NI(i)==NI(i-1))
		// 		min_max.push_back({(t(i)+t(i-1))/2,(NI(i)+NI(i-1))/2});
		// }

		// // combine values into row arrays
		// arma::Row<fltp> t_min_max(min_max.size());
		// arma::Row<fltp> NI_min_max(min_max.size()); arma::uword cnt = 0llu;
		// for(auto it=min_max.begin();it!=min_max.end();it++,cnt++){
		// 	t_min_max(cnt) = it->first; NI_min_max(cnt) = it->second;
		// }

		// // create mesh grid
		// const arma::Mat<std::complex<fltp> > z(
		// 	arma::repmat(arma::linspace<arma::Row<fltp> >(x1_,x2_,num_x_),num_y_,1),
		// 	arma::repmat(arma::linspace<arma::Col<fltp> >(y1_,y2_,num_y_),1,num_x_));

		// // calculate F
		// const arma::Mat<std::complex<fltp> > F = calc_F(z);

		// // create marchingcube algorithm
		// const cmn::ShMarchingCubesPr mc = cmn::MarchingCubes::create();
		// mc->setup_quad2isoline();

		// // std::cout<<NI_min_max<<std::endl;


		// // std::cout<<arma::real(F)<<std::endl;

		// // walk over local minima and maxima
		// for(arma::uword i=0;i<NI_min_max.n_elem;i++){
		// 	// get value 
		// 	const std::complex<fltp> z0 = create_aperture(t_min_max(i));
		// 	const fltp F0 = calc_F(z0).imag();

		// 	// create iso lines
		// 	const arma::field<arma::Mat<fltp> > Rl = cmn::MarchingCubes::combine_line_segments(
		// 		mc->polygonise(x1_,x2_,num_x_, y1_,y2_,num_y_, arma::imag(F), F0));

		// 	for(arma::uword i=0;i<Rl.n_elem;i++)
		// 		std::cout<<Rl(i).t()<<std::endl;
		// }

		
		// arma::real(F).eval().save("field.csv",arma::csv_ascii);

		return DFPolygon::create();
	}

	// bounding box
	arma::Mat<fltp>::fixed<2,2> DFMultipole::get_bounding() const{
		// allocate
		arma::Mat<fltp>::fixed<2,2> Mb;
		
		// bounding in x
		Mb.col(0) = arma::Col<fltp>::fixed<2>{x1_,x2_};
		
		// bounding in y
		Mb.col(1) = arma::Col<fltp>::fixed<2>{y1_,y2_};

		// return box
		return Mb;
	}

	// distance function
	arma::Col<fltp> DFMultipole::calc_distance(const arma::Mat<fltp> &p) const{
		return {0.0};
	}

	// validity check
	bool DFMultipole::is_valid(const bool enable_throws)const{
		return true;
	}


	// type string for serialization
	std::string DFMultipole::get_type(){
		return "rat::dm::dfmultipole";
	}

	// method for serialization into json
	void DFMultipole::serialize(Json::Value &js, cmn::SList &list) const{
		// parent
		DistFun::serialize(js,list);

		// store type
		js["type"] = get_type();

	}

	// method for deserialisation from json
	void DFMultipole::deserialize(
		const Json::Value &js, cmn::DSList &list, 
		const cmn::NodeFactoryMap &factory_list, 
		const boost::filesystem::path &pth){
		
		// parent
		DistFun::deserialize(js,list,factory_list,pth);

		
	}


}}